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Abstract. In this note, we provide simple convergence analysis for the algebraic multilevel itera- 
tion methods. We consider two examples of AMLI methods with different polynomial acceleration. 
The first one is based on shifted and scaled Chebyshev polynomial and the other on the polynomial 
of best approximation to on a finite interval with positive endpoints in the uniform norm. The 
construction of the latter polynomial is of interest by itself, and we have included a derivation of a 
3 term recurrence relation for computing this polynomial. We have also derived several inequalities 
related to the error of best approximation, which we applied in the AMLI analysis. 



1. Introduction 

The classical algebraic multigrid (AMG) algorithm was originally proposed in the early eighties 
[T3l [m [15] . From a practical point of view its description in [36] has made a great impact. In the last 
decade there has been a revival of research in the field of classical AMG [22], element-based AMG 
[181 [26l 129] . spectral AMG [21], AMG based on energy-minimizing interpolation and smoothed 
aggregation [38l dSl dHl [521 [5l], bootstrap AMG [16], adaptive smoothed aggregation [19] and 
adaptive algebraic multigrid [20]. While AMG methods can be viewed as matrix-based multigrid, 
there exist also various multilevel (ML) methods that utilize a sequenze of coarse-grid operators 
that arise form rediscretization of the continuous problem-based on its variational formulation-or 
from hierarchical basis transformation of the discrete problem. Multilevel methods evolved from 
two-level methods that have been introduced and analyzed for finite element elliptic systems in 
[H [11]. For sharp two-sided bounds on the convergence rate of two-level methods see [5^. A 
straightforward recursive extension of the two-level methods leads to the class of hierarchical basis 
(HB) methods for which the condition number of the preconditioned system in general grows 
exponentially with the number of levels, cf. [56]. That is why one often combines HB methods 
with certain stabilization techniques [50]. By using specially constructed matrix polynomials HB- 
ML methods can be made optimal. Known as (linear) algebraic multilevel iteration (AMLI) such 
methods have been introduced in the late eighties [T] [HI (371 [5T] , for nonlinear AMLI methods see 
also [5] [6l [32j [STJ [5l] . Recent works on algebraic multilevel preconditioning methods mainly focus 
on additive and parameter-free methods [3 [32l [l0| , aggregation-based preconditioners [121 [36l Il2] , 
and on extending their theoretic foundation and their applicability to various nonconforming and 
discontinuous Galerkin (DG) discretizations of highly ill-conditioned problems [21 [lOl HT] [231 (SI 
[251 [Ml [35]. 

Such problems arise for example in continuum mechanics when highly heterogeneous or porous 
media such as human bone tissue or geocomposites are considered, which have a complicated 
hierarchical structure with multiple characteristic length scales. In p3] preliminary results are 
shown on the construction of robust two-level methods for elliptic problems with extremely rough 
coefficients. 

The combination of AMLI and AMG techniques is a promising approach in view of increasing the 
robustness of multilevel solvers for this kind of applications. For instance it is known that standard 
AMG with piecewise constant interpolation yields an instable multilevel method in general even 
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though the related two- level method is stable, see, e.g., |51| . In such cases polynomial acceleration 
is a key tool for stabilizing the multilevel solution method. 

A theoretical comparison of algebraic multigrid and multilevel methods can be found in [51]. A 
common framework for their analysis has been established in form of the abstract theory of subspace 
correction methods |53^ [55] . A comprehensive exposition and analysis of various multilevel block 
factorization preconditioners, including AMG, has been presented in j51j . 

In this work we review the basic construction of multilevel methods as recursive extensions of 
two-level methods and consider polynomial acceleration techniques that can be employed in AMLI 
and AMG methods. In particular, we study the polynomial of best approximation to l/x on afinite 
interval [Amin , Amax], < Amin < Amax in the II • lloo (uniform) norm, derive recurrence relations and 
error estimates and show how these enter the convergence analysis and condition number bounds 
of the accelerated multilevel methods. 



In this section we consider a two-level iterative method for the solution of a system of linear 
algebraic equations 



where A £ ^ , is symmetric and positive definite, 5 G IR is a given right hand side. To 
describe a general two-level or multilevel multiplicative method, we denote V = IR^, and also 
introduce a coarse space Vh, Vh C V, Nh = dimVf/, Nh < N. In the following we will always 
assume that Vh = range(P), where P : IR"^^ i— ?• V and its matrix representation in the canonical 
basis of IR^^ is given by the coefficients in the expansion of the basis in Vh via the basis in V. 
Clearly P is a full rank operator, and its matrix representation is oftentimes called prolongation or 
interpolation matrix. 

2.1. Exact and inexact two-level methods and polynomial acceleration. We now consider 
a classical two- level preconditioning iteration, which utilizes a multiplicative preconditioner ~ 

Algorithm 2.1. Given x V the action B^^x is defined via the following three steps: 

1. Pre-smoothing: y = M~^a;. 

2. Coarse grid correction: z = y + PW^P^ {x — Ay). 

3. Post-smoothing: B~'^x:=z + M^^ {x — Az). 

To apply this algorithm, for example in the Preconditioned Conjugate Gradient (PCG) method, 
one needs to define an appropriate smoother M and a coarse grid preconditioner Bh- The coarse 
grid preconditioner is an approximation to the restriction of A on Vh- This restriction we denote 
by Ah and we have Ah = P^AP and BJ^^ ~ ^h^- To write down the closed form of B~^ we also 
need the symmetrization of M, denoted here with M. This symmetrization is defined in a usual 
manner, namely M satisfies 



2. Two-level and multilevel methods 



(2.1) 



Ax = b, 



{I - M-'^A) = {I- M-'^A) (I - M-^A) . 



A simple calculation then shows that 



(2.2) M = M{M + M'^ -A) M^, and M'^ = M'^ + M''^ - M' 



'^AM-\ 



Writing the three steps in the two-level algorithm in terms of x then leads to 



z 



M-^x + PBjj^P^ {x - AM'^x 



+ PBjj^P'^{I - AM~^ 
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and hence 

B-^x = z + A'r'^{x- Az) = Ar^x + {I -M-'^A)z 

= M~'^x + (/ - M^^A) \Ar^ + PBjj^P'^{I - AJVr 

= M~^x + (/ - A'r'^A)M-^x 



X 



+(/ - A)PBjj'P' {I - AM-')x 



X. 



= M"^ + (/ - A)PBh'P' {I - ANF 
Since this identity holds for all x we have that 

(2.3) = Kr^ + {I - M-'^A) PBJ^P'^ (I - AM-^) . 

The aim is to construct B such that the condition number B~^A is constant (independent of the 
size of A). We will always assume that the action of requires 0{N) operations. 
We distinguish two cases in choosing Bh'- 

(2.4) exact two level method if B^ = A^, where A^ = P^AP. 

(2.5) inexact two level method if Bjj ^ Ah- 

The exact two-level method is easier to analyze, in many cases has a uniformly bounded condition 
number, but is in general not practical, since Nh is usually taken proportional to N, and typical 
cases are Nh = N/2, Nh = N/4 or Nh = N/8. Then the size of the coarse level problem is 
comparable to the size of the original problem and hence the application of B~^ would require the 
solution of a coarse level problem which is of the same difficulty as the solution of (j2.ip . 

On the other hand, in a typical inexact method the action BJ^^ is defined recursively, by the same 
algorithm. Clearly, in such case, computing the action of B~^ on a vector requires 0{N) operations. 
Such method is known as multilevel V-cycle preconditioner. If in addition the condition number of 
B~^A is uniformly bounded, then we have an optimal method. However, depending on the choice 
of the coarse spaces Vh in the recursive application of the algorithm, one may also end up with a 
method for which this condition number depends on the size of the problem A^. Here coarse spaces 
is in plural, since recursive application of Algorithm 12.11 requires more than one Vh- A typical 
example for such behavior is when A results from finite element discretization of a second order 
elliptic equation with conforming linear elements and the coarse spaces Vh are chosen to interpolate 
piece- wise constant functions. In such cases one definitely wants to further improve the method by 
modifying Bh, without changing the coarse spaces. 

In our focus will be the Algebraic Multi-Level Iteration (AMLI) methods, which are constructed 
as follows: Suppose that we have an initial coarse level preconditioner Bh- This means basically 
that we have in hand a Bh-, and as an initial algorithm we consider Algorithm 1 2 . 1 1 with Bh instead 
of Bh- We would like to design now a Bh-, in such a way that the inexact method becomes close 
to the exact two-level method. In doing this, we would like to keep the number of operations in 
calculating the action of B~^ under control, and also to decrease the condition number of B~^A- 

In AMLI for a given initial Bh the coarse- level preconditioner Bh = BJ^^, is defined as 

(2.6) Bjj',^:=q,.^{Bjj'AH)Bjj\ 

where g^-i € Vu-i is an appropriately chosen polynomial of degree less than or equal to {v — 1) 
for some u >1. To give a motivation for doing this, we list here three examples, which may be not 
practical, but are instructive: 

• Set u = \ and qQ{x) = 1, which results in Bh = Bh, that is we have no change in the 
preconditioner. 
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Set u = Nfj and for a given non-singular operator Xh '■ Vh i— ^ Vh, define gNn-ii-^n) = ^h^- 
This results in the exact two- level method. Note that regardless of what form the initial 
coarse level preconditioner Bfj : Vh i— )• Vfj takes, as long as is invertible the resulting 
method will be the exact two-level method. 

Again, let u = Nh and in the notation of Algorithm l2.H th = P^{b — Ax). Let qNH-iiXn) 
be any polynomial such that 

qNH~l{XH)rH = Xjj^TH. 

As a result we get that BJj^^th = Aj^^rn, but the action of BJ^^^ depends on the argu- 
ment, and hence we have a nonlinear preconditioner. For a given th, such a polynomial 
)r/fcan be constructed by applying — 1) Preconditioned Conjugate Gradient it- 
erations to the solution of 

Ah^h = fH^ with preconditioner BJ^ . 

For V = Nh as we have in this example, we immediately obtain the exact two level method. 
Since applying (Nh — 1) PCG iterations may be expensive, one may take a fixed number 
of PCG iterations, namely u = 0(1) and the resulting method is known as nonlinear AMLI 
method. More details about nonlinear AMLI methods are found in 



2.2. Error propagation operators and polynomial acceleration. The two-level multiplica- 
tive preconditioner B~^ constructed in the previous section may also be used as an iterator. A 
two-level iterative method generates a sequence of iterates a;(j) (for a given initial guess £C(o)) 

(2.7) x^i_^^l)=x^i) + B~'^r(^i)=X(^i)+B-^{b-Ax^i)), i = 0,l,... 

If x^ is the exact solution to (12. ip and e^j) = — a;(j) is the error, then the error propagation 
operator E of the two-level iterative method (j2.7p is written as 

(2.8) e(,+i) =^e(,), E = I - B-^A. 

As we pointed out, an inexact two- level method is based on approximating the coarse grid 
correction for the exact two level method 

bh = Ajj^rn, with rn = P'^fi^i). 

Hence the error transfer operator for one coarse-level update Xfj^^ + P B~j^ P'^ r (j^^ is given by 

(2.9) Eh = I- PBJj^P^A 

where BJ^^ is the inverse of the initial preconditioner (the approximation to AJ^^) at the coarse 
level, which is used to build Bh- 

The following simple identity will be useful in deriving the error propagation relation between 
the coarse and the fine-grid level. 

Lemma 2.2. Let pi,{x) = l — qi,-i{x)x. Then the following identity holds for the error propagation 
operator E = I — B~^A, where B^^ is defined in Algorithm \2.1\ with B'^ = Bjj^^: 

(2.10) / - B~^A ={I - M-'^A) Pu{Eh) (/ - M'M) 
where 

(2.11) p^{xy.=p^{i - x). 

Proof. Let Qr be a polynomial of degree less or equal to r and Ah = P^AP denote the coarse-grid 
matrix. We will first prove that 

(2.12) Pqr{BjlAH)BjlP^A = qr{PBjl P^ A)PBjl P^ A. 
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It suffices to prove the identity (|2.12|) for the case qr{x) = x^, k < r. For k = the result is obvious. 
For k = 1 we obtain 

which is true since Ah = P^AP. Finally, assuming that ()2.12p holds true for k = j — 1 we have 

p [b-^AhY bjIp'^a = P {b-^AhY"^ BJj^P^APB-^P^A 

= {pbjj^p'^aY"^ pbj^p'^apbjYp'^a 

(2.13) = [PBjIp^AYpB'^^P'^A. 
Using (|2.3p we find 

I-B^^A = I - M~^A- {I - M-^A) PB^]^P'^ [I - AM"^) A 
= {I- M-^A) {I - M"^A) 
- (I - M-'^A) PBjI^P^A {I - M-^A) 

(2.14) = (I - M-^A) {l - PBjI^P'^A^ [I - M~^A) . 
Moreover, since Pv{x) = 1 — qu^i{x) x we can rewrite (j2.6p in the form 

Bh]u = {I-Pu{BjIAh))A-h^ 

= {I -{I- q,^i{BYYAH)) BjYAh) A~^ 

(2.15) = q,,^x[B-^ Ah)B~iI . 

Then by substituting (j2.15p in (j2.14p and using (j2.13p and finally (j2.9p we obtain the following 
representation of (12. Sp : 

l-B^^A = {I - M~^A){I - Pqu-i{B]j^AH)BJj^P^A){l - M-^A) 

= {I- M-'^A){I - qu-i{PBjlP^A)PBjlP'^A){l - M'^A) 
= {I- M-^A) (I - g,_i(/ - Eh){I - Eh)) (/ - M^^A) 
= {I- M-^A) {p^I - Eh)) {I - M-^A) 

(2.16) = {I- M~^A) PuiEn) {I - M^^A) 

□ 

We move on to describe the multilevel case and prove a condition number estimate, which 
depends on the approximation properties of Vh as well as on the estimates involving qi^-i{-). 

2.3. The algebraic multilevel iteration (AMLI) algorithm. We now focus on more imple- 
mentation and analysis details for multilevel methods, which are obtained by recursively applying 
Algorithm 12.11 In what follows we will denote by B^^^ a preconditioner for a finite element (stiff- 
ness) matrix corresponding to a A; times refined mesh {0 < k < i). We will also make use 
of the corresponding A;-th level hierarchical matrix A^''\ which is related to A^^^ via a two-level 
hierarchical basis (HB) transformation J^^\ i.e., 

(2.17) I^'^) = (jW)^A('^) JW. 

By aIj'^ and A\j\ 1 < i,j < 2, we denote the blocks of A^'^^ and A^''^ that correspond to the 
fine-coarse partitioning of degrees of freedom (DOF) where the DOF associated with the coarse 
mesh are numbered last. 
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The aim is to build a multilevel preconditioner B^^^ for the coefficient matrix ^(^^ := Ah at the 
level of the finest mesh that has a uniformly bounded (relative) condition number 



■^W-i^W) = 0(1) 



and an optimal computational complexity, that is, linear in the number of degrees of freedom N£ 
at the finest mesh (grid). In order to achieve this goal hierarchical basis methods can be combined 
with various types of stabilization techniques. 

One particular purely algebraic stabilization technique is the so-called Algebraic Multi-Level 
Iteration (AMLI) method, where a specially constructed matrix polynomial p^'^^ of degree Uf^ is 
employed at some or at all levels k = ko + 1,--- ,£. The AMLI methods have originally been 
introduced and studied in a multiplicative form, see Oil]. The presentation in this section follows 
Reference [32] . 

Starting at level (associated with the coarsest mesh) on which a complete LU factorization of 
the matrix A^^^ is performed, we define 

(2.18) 5^:=^^. 

Given the preconditioner at level k — 1 the preconditioner i?^'^^ at level k is defined by 



(2.19) 
where 

(2.20) 

and 
(2.21) 



_g(fe) i{k)j^(k)jj{k) 



^21 '-'1 





/ 



I 



A 



(k) 

12 



C 



(k) 



11 






^(k-1) 



Here we use the approximation 
(2.22) Z^^-^'^ := A^^^^U I 



(k)(r^(k~iy 



T{fc)^{fc)-l J(fe) 
^21 "-"11 



to the Schur complement S = A^''~^^ 
level stiffness matrix (stiffness matrix at level k 
hierarchical basis representation 



A\l' where ^('^-^^ 



A 



H 



^22^ is the coarse- 



1), which can be obtained from the two-level 



(2.23) 



^11 ^ 



12 

(k) Jik) 
22 



A'^ ' A 



(k) 
11 

w 

21 



T(fc) 

^12 



at level k, and p'-'^-* is a polynomial of degree satisfying 



(2.24) 



p(*^)(0) = 1. 



It is easily seen that (I2.22P is equivalent to 
(2.25) z(^-i)-^ = 

where the polynomial q^^^ is given by 



(2.26) 



x) 



1 - p^^\x) 
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We note that the multilevel preconditioner defined via (j2.19p is getting close to a two-level method 

when q^^\x) approximates well — in which case Z^^^^^ ~ ^(^^-i) ^ jj^ order to construct an 

X ^ 

efficient multilevel method the action of Z^^^^^ on an arbitrary vector should be much cheaper 

to compute (in terms of the number of arithmetic operations) than the action of A^^~^^ . Optimal 

order solution algorithms typically require that the arithmetic work for one application of Z^^~^^ 
is of the order 0{Nk~i) where N^-i denotes the number of unknowns at level k — \. 

A linear system with B^''\ an unknown vector v^^^ = {v^^ ^^^2^^ )^) ^ind right hand side vector 

d^^^ = {d^^^ , d'^^ )^ at level k can be written 



c['i> 



/ a 





(fe)- 



as 



-1 



A 



(fe) 



11 ^12 



V 



(k) 



Its solution involves the solution of a system 

at level (A: — 1) for some right hand side d^^~^\ Let 



(fc) 
1 

(fe) 







( uf^ \ 


\- 








then we have 














= I 


(2.27) 






= t 



11 

/ 



^ "-^11 ^12 



V 



V. 



(k) 

1 

(fc) 



Using (j2.25p we write (12.27P in the form 

(2.28) B^^-^^u'-^^ = gW(A('="i)i?(^-i)^ V^'"'^ 



where q^^\a 



+ a['''x + . . . + all^_^x''''-^ . Hence 



u 



(k) 



The following algorithm computes the solution of 

(2.29) i?Wt;W = dW. 

Algorithm 2.3. [Linear AMLI] (cf., [8]) 



8 



JOHANNES KRAUS, VERONIKA PILLWEIN, AND LUDMIL ZIKATANOV 



forward: 



(2.30) 



for A; = 1 to £ set ak-=0 
k:=i 

if o-fc = 1 

:= (jW)^dW 



w 



(fc) 



1 

4i ^1 



(fc-i) 



(2.31) 



(2.32) 



backward: 



rf(A:-l).=^(A:-l)^(fc-l) + 

end 

k:=k-l 

if A; > goto forward 

solve = for v^^^ 

k:=k+l 

if (7fe < Vk goto forward 



^12 ^2 



(2.33) 



(Tfc:=0 

if A; < ^ goto backward 



Here the vector u = {i'i,U2, ■ ■ ■ , ^i)^ defines the cycle, i.e., z^^ = 1 for 1 < k < £ corresponds to 
the V-cycle, and i^^ = 2 for 1 < k < £, I'l = 1, corresponds to the classical W-cycle. Higher-order 
stabilization or mixed cycles with varying polynomial degree are possible and sometimes preferable 
from a computational point of view. 

Algorithm ()2.3p is based on the multiplicative two-level preconditioner 



(2.34) 



C 
A 



(k) 
11 
(fc) 

21 



7ik) 

^12 



12 



' , and Z(^-i) for the inverse of the Schur complement 



It uses the approximations c[\^ for 
at level k, see (12:25]) . The action of the matrix polynomial qi^') (A^^'"^^ B^''^^^ ) in ([2:25]) on a vector 
is computed via Horner's rule, where the coefficients of q^^\x) are given by a^^^_^^, 1 < (t^ < i/^. 



The matrices A^2 -^21 off-diagonal blocks of the two-level hierarchical basis matrix 

A^*^^ at level k, and A^''~^^ is the matrix associated with the coarse grid (with respect to the 
coarse-grid nodal basis). 

2.4. Condition number estimates. Let us now study the spectral condition number ^(i?^^-' A*-^^) 
where B^^^ denotes the recursively defined linear AMLI preconditioner, cf. (j2.18p - (|2.22p . The pre- 
sented results are in the spirit of [31 [71 [8] and have the same recursive structure, that is, an estimate 
at level k involves the same type of estimate at level A: — 1. 

The basic assumption in the analysis of the multilevel preconditioner is an approximation prop- 
erty of the form 



(2.35) 



ei;'^v^A^'^v<v^C^''^v<ePv^A^''^v yv, k = l,2,...,i. 
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Additionally, it reasonable to assume that 

(2.36) 0<9o< eil'^ < ef^ < 01 < oo, l<k<l. 

The following lemma shows the basic induction step in the multilevel analysis. 
Lemma 2.4. // the multilevel preconditioner at level j = k — 1 satisfies the relation 

(2.37) pI^'^v^A^^^v < v^B^^^v < p'f^'^A^^^v 

then the preconditioner at level j = k satisfies the same relation with the constants 

max|l,rj^ | 

(2.38b) pf^ - 



mmjljTQ I 
where 

(2.39a) ri''"^^ = min xq^^\x), 

(2.39b) rf"^^ = max xq^''\x), 

and q^^^ can be any polynomial that satisfies q^^\x) > Vx G = [l/Pi'^ ^\ ^/ Pq ^^]- 

Proof. In order to prove the transition from level A: — 1 to level k we consider the matrix (j2.34p 
corresponding to the multiplicative two-level preconditioner at level fe, which can be factored as 



where 



cjf 
A(^-i) 



and the matrix L^^'^ is given by (j2.20|) . Now from (j2.25p we see that for any vector ^2 7^ and 
X2:={A^^~^'')~^/'^y2 we have 

yl{Ai'^-^))-^y2 

y^{^^^W^2 

a;^(A('^-l))V25(fe-l)-lg{fc)(^(fc-l)5(fe-l)-l)(^(fc-l))l/2^2 



X lXq^^\X)x2 



X^X2 

where X = ^A^k-i)^/^ ^{k-ir^ ^j^ik-i)^/^ ^^^^ ^^le relation 

(2.40) rt'\lZ^^-')y2 < ylA^^~%2 < rf^ yj Z(>^^%2 

holds with constants Tq'^ ^'^ and r^^^ defined in (j2.39p . Further, 

^Tc(k)^ _ yTD(^)y _ ylc[1^yi + y^A(''-^)y2 



(2.41) 



vTbWv yTDi^)y y^cf^y^ + yjz(^-i) 



y2 
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for all I' / and y = {yj, yj)"^ = L^^^v and thus, by using (12.40p . we find that for all v & V, 

(2.42) mm{l,r^Q~^'>}v^B^^'^v < v^C^^'^v < nia-K{l,rf'^^}v'^B^^^v, 
Finally, we combine the estimates 

(2.43) max , < max max , 
and 

(2.44) min — r-^min— ^ — ttt— < min — ^ — 77^, 

vy^o v^Cy'^'v v^O v^B^^'v v^OV^B^^'V 

with (|2.42p and (|2.35p and for all v £ V we obtain 

(2.45) 4x min{l, A^'^h B^'^^v < v^A^'^^v < 4t max{l, "'H v^B^'^^v, 

which proves the lemma. □ 

The following theorem provides an estimate of the (relative) condition number of the multilevel 
preconditioner b'^^\ 

Theorem 2.5. Assume that the approximation property \2. 35\l is satisfied on all levels k, 1 < k < i, 
which means that i2.36]) is true for some positive constants 6q and 61. Further, let < po < < 
pf^ < pi < 00 for all k = 1,2, ... Then the estimate 

(2.46) .{Bi^r^A') < ^^^l^W 
holds where 

(2.47a) ro = min xq'^''\x), 

a;e[pi~i,po"M 

(2.47b) ri = max xq^^\x). 

x6[pi"l,po~^] 

The condition number k{B^^^ A^) is uniformly bounded if 

(2.48) ^max{W < 

6'omm{l,ro} po 

Proof. Since 5^°) = implies Amin(S(°)~^A(°)) = Amax(S^°^"^^^°^) = 1 the bound (fOTD holds 
for J = and p^lf^ = p^^^ = 1. Then a repeated application of Lemma 12.41 finally results in the 
estimate (|2.46p for the multilevel preconditioner B'^^\ Clearly this estimate does not depend on 
the number of levels £ if (I2.48P holds. □ 

3. Best polynomial approximation to l/x in uniform norm 

In this section we prove a three term recurrence relation for the polynomials of best uniform 
approximation to l/x on a finite interval. We use symbolic computations to prove the recurrence 
relation. The specific form of these polynomials that we derive here is suitable for the AMLI 
analysis from the previous section and we provide such analysis in Section HI 

The polynomial of best approximation in uniform norm to l/x on a finite interval has been 
known for some time and can be found in classical texts on approximation theory, for example, 
see [Ml HSl d]. The problem of approximating x~^ has also been considered by Chebyshev [37] 
more than 100 years ago. A more recent work on this topic, which we will use in the derivation of 
the recurrence relation is 1281. 
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3.1. Notation. Our considerations are on a finite interval, [Amin > -^max ] 5 with < Amin < A^ax < 
00. Since our aim is to use the construction here for preconditioning of hnear systems, we introduce 

the condition number k = ^'^"^ > 1, and 



A. 



Ar 



+ A, 



Xr 



Ar 



Note that a > 1 and a > 0. The change of variables 

Am a x 



A. 



Ar 



+ Ar 



K+l 
K-1' 



2ax — a, 



maps the interval [Amax ; Amin] to [—1, 1]. The inverse map is 

x = — [t + a), and — = . 

2a X t + a 

We thus aim to find the polynomial of degree less than or equal to m of best approximation 

i'^ II ■ lloo f-i 11 norm of f{t) = , a > 1. We note that if Qmif) is the polynomial of best 

' ^ t + a 

approximation to + o) on [—1, 1], and the error of approximation is 

1 



■1,1] 



mm 



t + a 



Q 



then 



(3.1) Qmix) := 2aQmi2crx — a), and E = min 



oo,[-l,l] 



00,[A: 



2(7^[-l,l]. 



max 5/\min J 



are the polynomial of best approximation in || • ||oo norm on [Amin , Amax] and the error of approxi- 
mation respectively. 

We denote the Chebyshev polynomial of degree k by T^. For Tk{x) G Vk- we have 

1 \, 



We recall that 

and denote 
(3.2) 



^fe(*) = cos /c arccos(t), ts[— 1,1]. 



6 := a — \/g 



1 



T] 



Observe that 6 is the CG-convergence rate estimate, when Amax and Amin are the extreme eigen- 
values of a positive definite matrix, (see, e.g. M. R. Hestenes and E. Stiefel |27j). Note also that 

0<5 <1, = a+ \/a2 - 1, r? < and 5 = |r/|. 

With this notation in hand, we have the following identities, 



(3.3) 



2 t + a 2t — rj — rj ^ 



and directly from the expression for Tfc(^) given above, we also have 



(3.4) 



Next we give the best polynomial approximation to 



1 



t + a 



as it is given in [15] and [28] 
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Theorem 3.1. Let m > 1 be a fixed integer. The polynomial Qm G which furnishes the best 

approximation to in the L^o norm on [—1, 1] is 

t + a 

Qm{t) = — + — > V^Tj{x) - ziY^Tmix). 

T] — Tj ^ 7] — r] ^ [rj — T] ^j^ 

The error of best approximation is 

-Er„i 11 = mill 



Q 



oo,hl,l] " - 1 



t + a 

Proof See [45] or [28]. □ 

The following corollary is immediate and follows after a simple calculation. 

Corollary 3.2. -(/-E'm,[Aniin -Amax] error of approximation with polynomial of degree m on the 

interval [Amm,Amax], < Amin < Amax < oo then 

^m,[A.ni„,A„,ax] = 2<^'"~^^0,[yA— ,v/W]' 

where En, r a — n 1 is easily calculated to be 

I V '^min i V Amax J ^ 

1 



max J V \ / \ / \ 

^ VV'^min V^n 

For computations, the following result is useful, since it gives a straightforward way of computing 
the polynomial of best approximation. By means of algorithms for symbolic summation a simple 
closed form for the polynomial of best approximation stated above as well as a short recurrence 
relation can be found completely automatically. Nowadays there exist plenty of methods and 
implementations thereof that are up to complete these tasks [33]. It is in the nature of these 
algorithms that discovering new identities and relations means having a proof at hand at the same 
time. 

For our specific problem we choose to apply the Mathematica package SumCracker developed 
and implemented by M. Kauers |31j . We are going to use two of the features of this package: 
discovering (thus proving) closed form representations, and discovering (thus proving) recurrence 
relations for admissible sequences. Admissibility in this context means that the input for this algo- 
rithm consists of expressions that satisfy systems of difference equations, i.e., recurrence relations, 
with polynomial coefficients. This is certainly the case for Chebyshev polynomials. Informally 
speaking, the admissible sequences are translated into a difference ring context. In this difference 
ring computations are carried out using known techniques for commutative, multivariate polyno- 
mial rings, especially Grobner bases. The results of these computations can then be interpreted as 
statements about the original admissible sequences. For a detailed description of this and further 
algorithms contained in SumCracker we refer to \30\ [3T] . 

Theorem 3.3. The best approximation out ofVm to in the L^o norm on [—1, 1] is 

t + a 



Qm{t) — ^ , ( 1 — — — „_i\2 -^™+i(^) 



t + a\ (?? - If 
where 

Rm+l{t) = V~^T^+l{t) - 2Tm{t) + r]Tm-l{t). 

Furthermore for the polynomials of best approximation the following three term recurrence relation 
holds: 

(3.6) r]~^Qm+2{t)-2tQm+i{t) + vQrn{t) = -2, m = 0, 1, . . . 
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with 

Qo{t) = 

The error of approximation is: 



Qi{t) 



mm 



t + a 



Q 



oo,[-l,l] 



1 



Proof. The closed form can be obtained in an instance by an application of SumCracker's command 
"Crack" as follows: 

ln[l]:= Crack[SUM[r7''ChebyshevT[j,t],{j, 0,m - 1}]] 

-ry^+^ChebyshevTfm + l,t]+ t"'{2r]t ~ l)ChebyshevT[m, t\-T]t + \ 



0ut[l]= 



772 - 277* + 1 



From the sum representation given in Theorem 13.11 it is obvious that Qm{i) is indeed a poly- 
nomial of degree m. Certainly it is also an immediate consequence of the three term recurrence 
relation together with the initial values. This recurrence, however, can again be obtained completely 
automatic: 

ln[2]:= GetLinear Recurrence [p[m, t],In ^ m, Head — >■ p] 

Out[2]= p[m + 2] -277 - 77^p[77i] + 2trjp[m + 1] 

Here p[m, t] was defined as the polynomial of best approximation, the option "In" specifies the 
discrete variable along which we look for a recurrence relation and the option "Head" declares which 
variable shall be used for a concise output. The default head is set to "SUM". The initial values 
for the recurrence are easily computed as well, which completes the proof of the theorem. □ 



The next lemma gives an estimate on by a linear polynomial, which is used later to 

derive sufficient condition for the positivity of • 

Lemma 3.4. The following estimate holds for the polynomial Rm+i (t) defined in Theorem \3.3\: 
(3.7) -2{t + a) < R^+i{t) <2{t + a), t£[-l,l]. 

Proof. Recall that by the definition of rj and 6 (see (j3.2p ). we have that 77 < 0, and \rj\ = 6. Let 
t = cos a, for a G [0,7r]. Then we find that 

i2™+i(t) + 2t-r?-r?-i = i]-^{Tm+i{t) - 1) - 2{Tm{t) - t) + r]{Tm^i{t) + I) 

_i 2"^ + l m+l m — 1 

= — 2r? sin a + 4 sin asin a 

'2 2 2 

2 m — 1 
—27/8111 a 



(3.8) 
(3.9) 



_i I m+l 771—1 
-2ri sin a — 77 sin a 



^1/ m + l ^ 771,-1, 
26 M sin — - — a + (5 sin — - — a] > 0. 
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In an analogous fashion we obtain 

Rn,+i{t) - 2t + + T]-^ = i]'\Tm+i{t) + l)-2{Tm{t)+t) + r]{T^_i{t) + l) 

„i o'^+l m + 1 m — 1 

= 2r? cos a — 4 cos a cos a 

'2 2 2 

2 m — 1 
+2r/cos — - — a 

_i/m + l m — 1 

(3.10) = 2ri I cos — - — a — r/cos — - — a 



2 



2 



7Tt -\- 1 771 — 1 

(3.11) = -2(5"M cos— — a + Jcos— — a ) < 0. 



2 2 

Combining (j3.8p and (j3.10p and using 2t — — r/~^ = 2{t + a) yields the desired result. □ 

3.2. Algorithm for finding best polynomial approximation to 1/x. We now assume that 
we are given a matrix A, and bounds for its eigenvalues denoted by Amin and Amax • We will use 
the result from Theorem 13.31 to construct the polynomial qmi^) which approximates . We first 
write the recurrence relation for q^+iix) = 2aQm+i{2cTX — a): 

Qm.^i{2ax — a) = r][—2 + 2(2ctx — a)Qm{'^f^x — a) — r]Qm~i{2crx — a)]. 

Multiplying by 2a then gives 

qm+i{x) = r][-4a + 2a{2ax - a)Qm{2crx - a) - 2ariQ„T,^i{2ax - a)]. 

We then have the following algorithm in which the formulae are obtained by writing r/, a and a in 
terms of Amax and Amin • The polynomials qm{x) are written in terms of two parameters, x and 5, 
where 5 is defined in (j3.2p and for x we have 

^= ( A 1 /^ ^' 

(.V^max + V'^min ) 

Algorithm 3.5. Set /xq = 1/Amax and /ii = 1/Amin- 

1. Calculate the 0-th order polynomial go and the first order polynomial qi: 

9o(a;) = ^(Mo + w), and qi{x) = ]^{^/JI^ + y/JI{f - ^^^iix. 



2. For /c = 1, . . . , m — 1, g^+i written as a correction to q^ is computed as follows: 
Sfc+i(a;) = 4/io^i — ^ _ ^^^^^ _^ b\qi,{x) - gfc-i(x)] 



[l_g,(^)^]+52^^(^) 



gfc+i(x) = gfc(x) + Sfc+i(3;). 

The reason to write (?fc+i as a correction to is to show that such iterations look like iterations 
in a defect-correction method: First computing the residual [1 — qk{x) x], and then trying to correct 
it by adding an additional term. One can also easily see that for any initial go and gi, if the 
sequence q^ix) converges, then it converges to x~^ . In another word, choosing go and gi different 
from what they are above, will not generate the sequence of best approximations to x"^, but still 
this sequence will converge to x~^ . 

In the remaining sections we will show how the best polynomial approximation to 1/x can be 
utilized in multilevel methods, which results in efficient and often optimal order solution algorithms 
for symmetric positive definite systems of finite element equations, as they arise for instance from 
discretization of elliptic partial differential equations. 
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4. Application to multilevel methods 

Example 1. A classical choice of q^^^ is (I2.26P where p^^^ is a properly scaled and shifted Chebyshev 
polynomial satisfying the condition p^^\x) < Vt € [pi^^,po^^]- For degree two this polynomial 
is given by 

(4.1) p'^^\x)=p{x) = pQpx{^-s^{^-x 

which results in 

For this choice, 



.Po / VPi 
q^^\x) = q{x) = po + Pi - PoPia;- 



ro = min x q{x) = [x q{x)]\^^j_ = I, 

n = max xq[x) = \xq{x)\ L_ po+pi = — ■ 

Hence the uniform bound (|2.48p reads 

Oi {po + pif ^ pi 
6*0 4:PoPi ~ Po 
which requires a uniform approximation property such that 

W _ 



^0 - iPO+PlY + 

where we denoted hy R = —. 

Po 

Example 2. This example is analysis of AMLI method with the polynomial derived in Theorem 

13.31 Let Qm is the best approximation from Vm to in the Loo norm on [—1, 1]. We then set 

t + a 

qm{x) = 2aQm{2ax - a) 
where x = {t + a) /(2a). Recall again that from ()3.2p . we have rj < 0, and |r/| = 5. 

xqm{x) = 2ax Qm{2ax - a) = (t + a) Qm{t) 

^^■^^ -1 ^^"^ P m-1 ^(ziZ!^/? (A .Mill 

Combining (j3.7p with (j4.2p we obtain the estimate 



™ax^'e[pi-i,po J ' ^ 



maxie[_i,i] (l+ (^^pT)7(t + a)) 



ram^elpi-\po-']xqm{x) mintg[_i,i] (l - (^^f^(t + a: 

4(5'" 
(5-<5-i)2 



(4.3) < 



4<5" 



Tw(l + a) 



Substituting 6 = [^/k — l) / (\/k + l) and (1 + a) = {2k) /{k — 1) in the right hand side of ()4.3p we 
get an estimate that depends on the condition number k 

^4 4^ ™ax^g[pi-\po-i]a;gm(x) ^ 2 + ^""(At- 1) 

minxe[pi-i,po-i] 9n^(^) ~ ^ - (5'^(k - 1) 
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If the degree m of the polynomial is chosen large enough, that is, if 



2 + ^"^(k- 1) ^ 9o 
2-(5™(k- 1) - 6*1 -' 

\/s — 1 

for some number R > 1, where 5 = —= , or, equivalently, 

+ 1 

m > ^ — - 

log 5 

for given constants and 9i in (j2.35p then k is a uniform bound for the condition number k = 
i.e., K < K holds for any number of levels i. 
The best linear approximation to — on [pi~^ , oq^^] is given by 

X 

q{x) = ^(Vpo + ^TPif - PoPix. 

In this case the estimate (14. 4p shows that a uniform condition number bound R = — requires a 

Po 

uniform approximation property (j2.35p with 9q and 6i satisfying 

6*1 ^ k(1 + 2^^- R) 
Oq ~ 3 — 2^/r + R 

For comparison, if we would like to have k < 3 according to (j4.4p we need — < \/3, whereas the 

Chebyshev polynomial in Example [1] yields the same upper bound on the condition number already 

for— < — . 
^0 " 4 

Example 3. One interesting application of of the best polynomial approximation to — is in con- 

X 

structing a smoother for any AMG (or AMLI) method. In order to use q'm(^) as smoother for A, 
we have to make sure that qm{x) > for all x in (0, A] where A is an upper bound on the spectrum 
of A. In practice, when estimates on Amin (A) and Amax (A) are not available one can use 



(4.5) A=P| 



and construct the best approximation qmix) to — on the interval [— , A] where /i > 1 is a constant. 

X _/i 

Note that qm{x) > for all x in (0, A] if qm(x) x > for all x £ [— , Al. The latter condition holds 
true if 

(4.6) 1- ^l_'^_,^^ Rm+i{t)>0. 
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We first estimate below the left side of this inequality as follows 
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> 1 
= 1 

> 1 





— ^j^ 








-(5-1)2 




4(5"^ 




-(5-1)2 








-5-1)2 



a - 1 



where we have used that 5 = a — \/a?'^-l and (5 — 5 ^)^ = 4(a^ — !)■ As a consequence, a sufficient 
condition for the inequality (|4.6p to hold then is given by 

(4.7) (^Y< ' 



An upper bound on the damping factor of the smoother on the error components that correspond 
to eigenvalues in the interval [— , A] is as follows 

(4.8) p^l-q^^A)A)<^^ ^^"^ 



For example, if /i = 4 the condition ()4.7p holds for all m > 1. Choosing m = 2 in this case 

results in a damping factor of at most — for all error components in the span of eigenvectors whose 

6 

eigenvalues are in the interval [^)^]- 

If ;U = 8 the condition (j4.7p for positivity of qm{x) is satisfied for all m > 2. Choosing m = 3 in 

this case the estimate of the damping factor (right hand side of ()4.8p ) is given by 2 ( — ^ 1 ~ 

0.381276 for the interval [■^)^]- To calculate the damping factor in other cases is also rather 
straightforward. 



5. Correspondence and differences between the standard AMLI methods and AMG 

Finally, we would like to comment on the correspondence and differences between the "smoother" 
M and the prolongation used in AMG and the ones used in the standard AMLI method. The 
standard AMLI method, as introduced in [TJ E] uses f-smoothing only, i.e., 

'^Cii — CiiAuC^^ 




(5.1) 
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One possible implementation (that corresponds to the multiplicative variant) yields the precondi- 
tioner 



B- 



-1 



+ 






/ 



Z-' [-A2iC^,\l] 





+ 

/ 




Cil^n -Cf/Ai2 
II 
(5.2) = M"^ + M-'^A) PZ-^P^ (/ - AM-^) , 

which uses the symmetrized smoother (jS.ip . Then, in order to write B^^ in the form ()2.3p one has 
to choose 





/-^iiCf/ 



(5.3) 



P 



We notice that the simple form of (|5.3p is related to the fact that the AMLI preconditioner is 



defined for the hierarchical two-level matrix A = j'^AJ, which contains already the coarse-level 
matrix sub-matrix in its lower right block, i.e.. 

Ah = [0,1] A 




/ 



This, however, is in accordance with the Galerkin relation Ah = P^AP as is used in AMG methods. 

We conclude that the polynomial acceleration techniques that we studied here can be exploited in 
various ways in the implementation of AMLI and/or AMG based preconditioners, since both can be 
viewed as inexact two- level methods as we described in Section 12.11 The performance of both type 
of methods crucially depends on the particular choice of the polynomial q^-i^-) in equation (j2.15p . 
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